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Abstract 

Autopsy studies of adults dying of non-cancer causes have shown that virtually all of us possess occult, cancerous lesions. 
This suggests that, for most individuals, cancer will become dormant and not progress, while only in some will it become 
symptomatic disease. Meanwhile, it was recently shown in animal models that a tumor can produce both stimulators and 
inhibitors of its own blood supply. To explain the autopsy findings in light of the preclinical research data, we propose a 
mathematical model of cancer development at the organism scale describing a growing population of metastases, which, 
together with the primary tumor, can exert a progressively greater level of systemic angiogenesis-inhibitory influence that 
eventually overcomes local angiogenesis stimulation to suppress the growth of all lesions. As a departure from modeling 
efforts to date, we look not just at signaling from and effects on the primary tumor, but integrate over this increasingly 
negative global signaling from all sources to track the development of total tumor burden. This in silico study of the 
dynamics of the tumor/metastasis system identifies ranges of parameter values where mutual angio-inhibitory interactions 
within a population of tumor lesions could yield global dormancy, i.e., an organism-level homeostatic steady state in total 
tumor burden. Given that mortality arises most often from metastatic disease rather than growth of the primary per se, this 
finding may have important therapeutic implications. 
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Introduction 

Almost all of us carry small tumor lesions that for many will not 
progress to symptomatic disease. Indeed, as evidenced in autopsy 
studies for adults without pre-established cancer such as [1,2], 
occult lesions are present in most healthy adults. Nielsen et al. [3] 
found that, out of 1 1 0 women cases, among which only one had 
been previously treated for breast cancer, 22% had at least one 
malignant lesion. Moreover, 45% of these had multicentric lesions. 
Similar results have been reported for prostate cancer in men [4] . 
For thyroid cancer, autopsy results [2] showed a prevalence rate of 
99.9% for occult carcinomas, while incidence of thyroid cancer is 
only 0.1% [5]. 

To explain these results, it is necessary to understand the tumor 
dormancy phenomenon. Tumor dormancy [6,7] is defined by 
stable or very slow tumor growth. It can happen at the cellular 
level as a malignant cell remaining quiescent for a long period 
before awakening, but here we focus on the mm-scale lesions such 
as have surfaced in the several remarkable autopsy studies 
discussed, i.e., tissue-level tumor dormancy. Although the sizes 
of these dormant tumors remain almost constant, it is not due to a 
cessation in cell proliferation, but rather to increased apoptosis 
that leads to a near zero net growth rate [6-8] . Clinically, tumor 
dormancy has been observed in breast cancer [3,9-1 1], melanoma 
[12] and prostate cancer [4], among many others [6]. Dormancy is 
particularly relevant to the situation where secondary tumors 
(metastases) remain small and undetectable for extended periods. 



Various explanations have been proposed for tumor dormancy, 
among these being the achievement of a balance between 
stimulation and inhibition of angiogenesis [7,13,14]. This mech- 
anism offers one explanation for how secondary tumors may be 
suppressed to a near-dormant state by the primary; a phenomenon 
known as 'concomitant resistance' [15,16]. In fact, a number of 
explanations for the concomitant resistance phenomenon have 
been suggested, as well summarized by Chiarella et al. [16]: 1) 
monopolization of certain resources by the primary tumor that 
deprives secondary tumors of materials needed for growth, 2) 
primary tumor-induced enhancement of immune suppression of 
small secondary tumors (concomitant immunity), .3) anti-prolifer- 
ative molecules released by the primary tumors and 4) release of 
angiogenesis inhibitors by the primary tumor into the blood 
circulation resulting in inhibition of vascular development at 
secondary sites. Nevertheless, although a distant impairment of 
metastatic growth by a primary tumor has been recognized for 
over a hundred years [17], and has meanwhile been informed by 
various preclinical [18-22] and clinical [9,23-25] studies, it 
remains poorly understood. 

However, because of evidence that concomitant resistance 
happens in immune-deficient mice [21] and considering the large 
and unequivocal body of support for the role angiogenesis 
inhibition plays in the maintenance of tumor dormancy [8,26- 
30] and the "angiogenic switch" [31] in escape from dormancy, 
our focus here will be on the last theory. Angiogenesis, the process 
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of creating new blood vessels and developing a supporting vascular 
network, was shown by Folkman [32] to be critical for tumor 
growth. Indeed, without development of new blood vessels, a 
malignant neoplasm cannot grow further than about 2 to 3 mm in 
diameter, due to nutrient supply limitations [32]. This process is 
regulated by the release from cancer cells of stimulatory growth 
factors, such as vascular endothelial growth factor (VEGF), that 
induce proliferation, migration and maturation of surrounding 
endothelial cells, as well as the production of angiogenesis 
inhibitory factors that act to curtail endothelial expansion [33]. 
As an example, in 1994 when examining the growth of Lewis lung 
carcinoma in a syngeneic murine tumor model, O'Reilly et al. [26] 
discovered an endogenous molecule having an inhibitory effect on 
angiogenesis, which they called 'angiostatin', followed soon by the 
discovery of 'endostatin' [27]. Endogenous anti-angiogenic mol- 
ecules were also evidenced in human cancer, an example being 
thrombospondin- 1 [29]. Overlaying the ability of tumors to 
stimulate vasculature, the discovery of their ability to also inhibit it 
[34] allows for the possibility that tumors may indirecdy control 
their own growth [14,33,34], perhaps as a vestige of normal organ 
growth control. Further, inherent to this self-control notion, if the 
inhibitors were longer-lived and thus more persistent in the 
circulation, they could have the collateral effect of suppressing 
angiogenesis and growth at distant metastatic sites as the tumor 
mass gets large [14]. Indeed, the half-life of angiogenesis 
stimulators has been reported to be on the order of minutes for 
VEGF [35], while that for angiogenesis inhibitors is on the order 
of hours [26,29]. 

Amidst these developments, there have been a number of efforts 
to take numerous complex mechanisms of cancer biology into 
account in mathematical models (see [36] for a review), but very 
few of these models have had the aim of describing metastatic 
development, despite metastasis being the main cause of death 
from cancer [37]. Indeed, while the cure rate of cancer before 
appearance of metastases is about 90% for all cancers combined, it 
falls to just 15% when distant metastases are present at diagnosis 
[16]. As far as we know, modeling efforts in this direction can only 
be found in the work of Liotta and coworkers [38], and more 
recently in a few stochastic models [39-42] describing progression 
through the different stages of the metastatic process (cell 
detachment, intravasation, survival in the blood, extravasation, 
settling in a new environment), and in one notable dynamic model 
[43]. In this last case, Iwata and coworkers [43] proposed a 
quantitative formalism for the development of metastatic colonies, 
which was of great potential interest as it is was designed to 
describe the size distribution of the metastases, allowing thus to 
distinguish between micro-metastases and larger lesions. However, 
this model does not take angiogenesis into account. We therefore 
decided to theoretically combine this work with that of Hahnfeldt 
et al. [14] for tumor development under angiogenic control, along 
with a mathematical model developed for the growth and 
dissemination of a metastatic population [44,45]. The goal we 
realized was a new global formalism that integrates local 
stimulation with systemic inhibition of angiogenesis by a circulat- 
ing factor produced by each lesion in a population of tumors, to 
provide insight into the development of the entire tumor/ 
metastasis system. 

Methods and Results 

In silico model - derivation and implementation 

The global philosophy of the model we propose is to consider 
the development of cancer disease at the organism scale, by 
describing the colonization and dissemination of a population of 



secondary tumors (metastases), in parallel with the growth of the 
primary lesion, taking into account organism-scale signaling 
interactions amongst these various tumor sites. The impetus for 
this viewpoint comes from Iwata et al. [43] where the authors 
derived a structured population model for describing the 
metastatic colonies represented by a density structured in size 
(volume). This model consists of a linear transport partial 
differential equation with a nonlocal boundary condition of 
renewal type. It has been further mathematically studied in 
[46,47], in particular to develop efficient numerical methods for 
discretizing the problem. 

A major limitation of this approach, though, is that it does not 
take angiogenesis into account, although this is a fundamental 
process of tumor development that cannot be neglected, partic- 
ularly if we want to study the effects of clinical angiogenesis 
inhibition. However, by combining the approach of Iwata et al. 
[43], arguably the first dynamical model for metastatic develop- 
ment, along with the model of Hahnfeldt et al. [14], which is the 
first to consider angiogenic homeostatic control of tumor growth, 
we developed in previous work a hybrid construct that integrates 
the angiogenic process into the growth of each tumor [44,45,48]. 
Since this model was written at the level of the organism, it was 
considered a suitable framework to adapt to the problem of 
analytically describing the consequences of systemic inhibition of 
angiogenesis (SIA). 

The result is a model for tumor growth control that takes into 
account the local and systemic actions of angiogenesis regulators. 
It integrates the ability of tumor lesions to locally stimulate 
angiogenesis while simultaneously inhibiting angiogenesis globally, 
and is fitted to preclinical data. Information on the behavior of 
metastases is inferred from the estimated parameters. Simulations 
of the cancer history are performed, which provide a detailed 
description of the distribution of predicted metastatic lesion sizes. 
The biological hypothesis of a global dormancy state of self- 
inhibiting tumors is then tested, and corresponding ranges of the 
inhibitor production rate identified. 

A schematic view of the new formalism we propose is presented 
in Figure 1. The main feature added to the previous model 
(Benzekry [48]) is a new variable representing the circulating 
concentration of an endogenous angiogenesis inhibitor, standing in 
for all possible inhibitory molecules (examples being endostatin, 
angiostatin or thrombospondin- 1) impacting on the growth of each 
tumor. As a general modeling principle, we sought to be 
parsimonious and describe the major dynamics of the system with 
as few parameters as possible to assure each dynamic introduced 
carries its proper burden to explain the data. 

Mathematics of tumor growth and systemic inhibition of 
angiogenesis 

Our construct considers primary tumors and their metastases to 
be distinct lesions whose states are described by two traits: volume 
V and carrying capacity K. The primary tumor state is denoted 
(V p (t), K p (tf). The model's main variable is p(t, V, Kj, the 
physiologically structured density of metastases having volume V 
and carrying capacity K at time /. The term density means that the 
metastases are assumed to exist in a continuum of sizes and 
carrying capacities and that the number of tumors between 
volumes V\ and V<i and carrying capacities between K\ and K 2 is 
given by p(t, V,K)dVdK. We assume that the dynamics of 

each tumor's state are governed by a growth rate for (V, Kj, 
denoted by the vector G( V, K; Vp, p), that dissemination of new 
metastases is driven by a volume-dependent emission rate /?(Pj, 
and that the repartition of metastases at birth is given by jV(F, Kj. 
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The precise expressions of these functions will be described below. 
We consider some fixed final time 7" and a physiological domain Q 
for the possible values of (V, Kj, defined as Cl = (V 0 , +°o)x(0, +00) 
where the distribution of metastases has its support, which means 
that metastases have size bigger than the size of one cell V 0 and 
non-negative carrying capacity. In the formula below, the vector 
v(V, Kj stands for the external unit normal to the boundary d£l of 
the domain Q. The notation 3Q + stands for the subset of the 
boundary where the flux is pointing inward, i.e. where G(V, K; Vp, 
p) • \{V, Kj<0. The map (V, K)\->p°{V, K) denotes the initial 
distribution of the metastatic colonies. 

Overall, the model we arrived at is a nonlinear transport partial 
differential equation of renewal type with a nonlocal boundary 
condition. 



d t p + div(pG) = 0 (0,T)xQ 
- G( V,K; V p ,p)»v( V,K)p(t, V,K) = N( V,K) 



> (1) 



/?( V)p(t, V,K)dVdK + P( V p (t)) \ (0, T) x dQ + 

! J 

p(0,V,K) = p°(V,K) Q 



We now make precise the expressions of the various coefficients 
of the model; in particular how the growth rate G is affected by the 
total population of tumors represented by p. We assume that all 
the tumors (primary and secondaries) share the same growth 
model but have different parameters, due to the different sites 
where they are located. However, within the population of 
metastases, all tumors are assumed to grow with the same 
parameters. The growth velocity of each tumor is given by a vector 
field G(V, K; Vp, p). Following the approach of [14] we assume 



G(V, K-V p ,p)-- 



K 



aV\n 



,Stim(V, K)-Inhib(V, K; V p ,p) 



In the previous expression, the first line is the rate of change of 
the tumor volume V (where a is a constant parameter driving the 
proliferation kinetics of the cancer tissue) and the second line is the 
rate of change of the carrying capacity K. The main idea of this 
tumor growth model is to start from a gompertzian growth of the 
tumor volume (or any carrying capacity-like growth model [49]) 
and to assume that the carrying capacity A" is a dynamical variable 
representing the tumor environment limitations (here limited to 
the vascular support) changing over time. The balance between a 
stimulation term Stim(V, Kj and an inhibition term Inhib(V, K; Vp(t), 
p(t, V, Kjj governs the dynamics of the carrying capacity. For the 
stimulation term we follow [14] and assume 

Stim(V, K) = bV, 

where the parameter b is related to the concentration of angiogenic 
stimulating factors such as VEGF. This last quantity was derived 
to be constant in [14] from the consideration of very fast clearance 
of angiogenic stimulators [35] . 

For the inhibition term, Hahnfeldt et al. [14] only considered a 
local inhibition coming from the tumor itself. Our main modeling 
novelty is to consider in addition a global inhibition coming from 
the release in the circulation of angiogenic inhibitors by the total 
(primary + secondary) population of tumors. The following is an 
extension of the biophysical analysis performed in [14]. Let us 
consider a spherical tumor of radius R inside the host body. The 
host is represented, for simplicity, by a single compartment of 
volume Vj in which concentrations are assumed spatially uniform. 
Let n(r) be the inhibitor concentration inside the tumor at radial 
distance r. Let the intra-tumor clearance of inhibitors, known to be 
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slow, be (approximately) zero [14]. At quasi steady state, n(r) then 
solves the following diffusion equation: 



n"(r) 



2n'(r) 



+ 



P_ 
D 2 



= 0. 



where p is the inhibitor production rate and D 2 is the inhibitor 
diffusion constant. This equation has the boundary condition 
n(R) = i(t; Vp, p(t, V, Kj), where the expression on the right 
represents the systemic concentration of the inhibitor resulting 
from a primary tumor volume V p and secondary tumors of density 
p at time t. Solving this equation (using that n(0)<+°°) we obtain 



dt V K, 



G p (V p ,K p ; V p ,p{t,V,Kj) 



(4) 



where G p has the same expression as G, except that the parameters 
a p and bp (the values of a and b that are associated with the function 
G for the primary) may be different from a and b associated with 
metastases, in those cases where the primary and metastases are 
presumed to have different growth kinetics. The inhibitor 
production rate p and effect of the inhibitor e are assumed to be 
the same for the primary and secondary tumors, which implies 
same value also for d in view of formula (3). 



n(r) 



i + 



P 
6D 2 



(R 2 



From this expression we compute the mean inhibitor concen- 
tration in the tumor to obtain 



Inhib{V,K; V p ,p) = e(i + j^R 2 )^ = 



where e is a sensitivity coefficient. For i(t; Vp, p(t, V, Kj), considering 
that the total flux of inhibitors produced by a tumor with volume V 
is /> 17 and assuming that the inhibitor production rate is the same in 
all the tumors, we have 



di 

V d - = P V p + 



\pVp(t,V, 



K)dVdK - kV d i, 



where k is an elimination constant from the blood circulation. 
Setting I(t; Vp, p(t, V, Kj) = Vjfa V p , p{t, V, Kj), we get 



dl_ 
dt 



= PV P + 



\pVp{UV, 



K)dVdK - kl, 



which has an initial condition that in significant cases may be set to 
zero, i.e., I(t=Q) = Q. Overall, the explicit expression of the 
metastases growth rate is given by 



G(V,K; V,,p) = 



hV-dVTK-elK , 



(2) 



where e = — and 
Vd 



d = eV^flY 



\5D 2 \4nJ 



(3) 



Note that we retrieve here the local term dV from the analysis 
of [14]. Our analysis results in an additional global term el that 
captures the effect of systemic inhibition of angiogenesis. 

For the primary tumor, we assume the same structural growth 
model. The dynamics of (V p , K p ) are thus given by 



Metastatic dissemination 

There is no clear consensus in the literature about metastases 
being able to metastasize or not [50-52]. However, we argue here 
that cancer cells that have acquired the ability to metastasize 
should conserve it when establishing a new site. Moreover, since 
metastases remain undetectable for an extended time [50-52] (in 
particular because tumors could remain dormant for some time), 
the absence of clear proof in favor of metastases from metastases 
could be due to the short duration of the experiments compared to 
the time required for a second generation of tumors to reach a 
visible size. Here we are interested in long-time behaviors and, 
although metastases from metastases could be neglected to a first 
approximation, we think this second-order term is relevant in our 
setting and chose to include it in our modeling, in light of some 
clinical evidence supporting second-generation metastases [53]. 

Successful metastatic seeding results when one malignant cell is 
able to overcome various adverse events including: detachment 
from the tumor, intravasation, survival in the blood/lymphatic 
circulation, escape from immune surveillance, extravasation, 
survival in a new environment (see [54] for more details about 
the biology of the metastatic cascade). Here, we regroup all these 
events into one emission rate ji(V, K), quantifying the number of 
successfully newly created metastases per unit of time. We assume 
very small metastases do not metastasize because they do not have 
access to the blood circulation, accounted for here by including a 
threshold V m below which tumors do not spread new individuals. 
V m is taken here to be 1 mm 5 as an approximation of the volume 
at which the angiogenic switch happens [32]. Apart from the 
addition of this threshold, the expression chosen for /? is the same 
as that used by Iwata et al. [43]: 



-{ 



mV a 
0 



otherwise 



(5) 



where m and a are coefficients quantifying the overall metastatic 
aggressiveness of the cancer disease. The parameter m represents 
the intrinsic metastatic potential of the cancer cells, and a 
represents the microenvironmental component of metastatic 
dissemination. It lies between 0 and 1 and is the third of the 
fractal dimension of the tumor vasculature, assumed here to be the 
same for all tumors. For instance, if vasculature develops 
superficially, then a = 2/3, whereas for a fully penetrating 
vasculature, the value would be a=l. We here assume the 
dissemination rate depends only on the volume because simula- 
tions revealed that adding a monotone dependence on A" did not 
improve the flexibility of the model even while adding at least one 
parameter, contrary to the parsimony principle. 

Stating a balance law for the number of metastases when they 
are growing in size gives the first equation of (1). The boundary 
condition, i.e. the second equation of (1), states that the entering 
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flux of tumors equals the newly disseminated ones. These result 
from two sources: spreading from the primary tumor, modeled by 
the term flVp(t), and second-generation tumors coming from the 
metastases themselves, described by the term 
J P(V)p(t,V, K)dVdK. The map (V, Kj^MV, K), where (V, K) e 
a 

dQ., stands for the volume- and carrying-capacity-dependent 
distribution of metastases at birth. Assuming that newly created 
tumors all have the size of 1 cell, denoted by V$, and some initial 
carrying capacity, denoted by K 0 , we have 

N(V, K) = V,Jt)=(>W ( 6 ) 

i.e. the Dirac distribution centered in (V 0 , K 0 ). We have previously 
discussed how this form can be deduced by passing to the limit 
from an absolutely continuous density [55]. An important feature 
of this model, in contrast to previous no-SIA models, is that we 
allow metastases to exit the domain by imposing the boundary 
condition only where the flux points inward and letting tumors exit 
the domain in the opposite situation. In view of expression (2), this 
occurs when the carrying capacity Kis less than the volume of one 
cell, V Q , i.e. when global inhibition is strong enough so that tumors 
can cross the line K= V Q , which is the case when G(Vo, 

2 

Fo)*(0, 1) = b Vo — dVtf> — elVo < 0. These tumors are then removed 
from the population, corresponding to death caused by nutrient 
deprivation. 

From the solution p of the model (1,3-6), biologically relevant 
macroscopic quantities can be defined, such as the total number of 
metastases N(t)= [ p(t, V, K)dVdK, the total metastatic burden 
h 

M(t) = [ Vp(t, V, K)dVdK, or the mean size of the metastases 
h 

M{t) 

my 

Solution-finding 

To approximate the solutions of the problem (1,3-6) we adapted 
a numerical procedure previously developed for the model without 
SIA in [45,55]. It is a Lagrangian scheme based on the 
straightening of the characteristics of the transport equation. We 
then used an Euler method for discretization of the characteristics 
and computation of the primary tumor ordinary differential 
equation. The integral in the boundary condition was computed 
using the trapezoid approximation method. 

Parameters surmised from existing preclinical data 

Data on metastatic development are not common in the 
literature, especially for micrometastases or dormant tumors, since 
these measurements are technically difficult to obtain. Even more 
difficult to find are data quantifying systemic inhibition of 
angiogenesis. For our purpose we use data from Huang et al. [56] 
that do not explicitly deal with systemic inhibition of angiogenesis 
nor global dormancy, but where number and mean size of 
metastases at the end time (T = 32 days) are available, together 
with primary tumor growth kinetics. The cell line used in this work 
is a spontaneous mouse breast cancer line 4T1, known to be highly 
metastatic with relatively slow primary tumor growth. Cells were 
injected subcutaneously (10 5 cells) in BALB/ c mice. As shown in the 
following, our model was able to explain these experimental data. 

Values of the parameters were fixed either by direct extraction 
from the literature, heuristic derivation, or by fitting the model to 
the data from [56]. For the preclinical data that we used, 
metastases actually develop to symptomatic volumes and did not 
evidence manifest global inhibition. Hence for fitting of the model, 



we consider SIA as being negligible and take /= 0. The parameter 
estimation we performed here is only intended for estimation of 
growth and metastatic spreading parameters. The assumption of 
negligible SIA was found a posteriori to be adequate for description 
of the data from [56], because adding an SIA term with reasonable 
parameter values did not have any impact on the model 
simulations, in the framework of the experiment from [56]. In 
the context of no SIA, there is no impact of the metastases on the 
primary tumor and we could separately fit the primary tumor 
growth and the metastatic development. This approach (compared 
to a global fitting of all the parameters together) further reduces 
the parameterization of the model and allows for more stable and 
biologically relevant parameter estimation. Indeed, only two 
degrees of freedom were used to fit the primary tumor growth 
(seven time points) and two for the data on metastases (two 
measurements). The meaning, units and values of the model 
parameters resulting from the whole estimation procedure are 
summarized in Table 1. 

In Hahnfeldt et al. [14], values for the elimination rates k and 
efficacy constants e for two endogenous angiogenesis inhibitors, 
endostatin and angiostatin, were estimated by fitting tumor growth 
data of mice that received injections of these anti-angiogenic 
agents. We focus here on angiostatin, and use values for the agent 
efficacy e and the elimination rate k from the blood circulation 
reported in Hahnfeldt et al. [14], applying these to a 20 g mouse. 
This value gives a half life for angiostatin of 1 .8 days, which is 
consistent with the value of 2.5 days that can be found in the 
literature [26]. 

O'Reilly et al. [26] showed that injection of 12.5 |U.g per day of 
recombinant human angiostatin reproduces the systemic inhibi- 
tion due to a primary tumor removed when it reached the size of 
1500 mm 3 . An approximation of the production rate in their 
12 5 

setting is -P* J^qq x 10~ 3 «8.3 x 10~ 6 mgmm~ 3 day _1 . For the 

value of V d we argue that the blood volume of a mouse is about 
1.2-1.6 cm' 5 per 20 g body weight. Taking an approximate value 
of 1.4 cm 3 and assuming that the interstitial (extracellular) space 
fills 30% of the extravascular space (in agreement with measure- 
ments of the fraction of volume occupied by cells), summing the 
interstitial space and blood volume gives 6.98 cm 3 . Hence we took 
Vj = 7000 mm' to be an approximation of the distribution volume. 
For the diffusion coefficient of angiostatin, D , we used a value 
1.56 mm 2 day 1 taken from the literature [57]. Based on these 
values and the formula (3) for d derived in the modeling section, we 
were able to heuristically compute an approximation of the 
parameter d as </~0.0717 mm~ 2 day _1 . In the following, we fixed 
d and dp to this value, which allowed us to reduce possible 
indeterminacy in the parameter estimation for the growth model. 

When reproducing the experiment of Huang et al. [56], we 
fixed the initial size of the primary tumor to be Vp(0) = 0.1 mm 3 
(corresponding to 10 5 cells, i.e. the number of cells injected into 
the mouse) and arbitrarily set the initial carrying capacity of the 
primary tumor to K p (0) = 200 mm 3 . Metastases were assumed to 
start with initial size Vq = 1 cell = 1 0 6 mm 3 and initial carrying 
capacity K 0 = 1 mm 3 (assumed to be an approximation of the 
maximum reachable size without angiogenesis [32]). For meta- 
static emission, we considered a superficial vascular development 
and took a = 2/3, following what was estimated in Iwata et al. [43] 
from clinical data. 

Fits to the data 

Parameters cip and bp were obtained by minimizing the sum of 
squared errors between the tumor growth model simulation and 
the primary tumor growth data from [56]. Least squares 



PLOS ONE | www.plosone.org 



5 



January 2014 | Volume 9 | Issue 1 | e84249 



Dormant Metastases from Angiogenesis Inhibition 



Table 1. Values, units and meaning of the model parameters. 



Parameter 


Value 


Unit 


Meaning 


Rationale 


a P 


0.154 


day 


PT cells proliferation 


Fit PT 


b„ 


16.7 


day 


PT angiogenic stimulation 


Fit PT 


d P 


0.0717 


—2 j —1 

mm day 


PT angiogenic local inhibition 


H 


a 


0.154 


day 


Met cells proliferation 


Fit Met 


b 


12.5 


day 


Met angiogenic stimulation 


Fit Met 


d 


0.0717 


mm day 


Met angiogenic local inhibition 
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PT = Primary Tumor. Met = metastases. I = global amount of angiogenic inhibitor in the blood. H = heuristic derivation. 
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minimization was performed using the trust region reflective 
algorithm implemented in Matlab (Matlab 2009b, The Mathworks 
Inc.). We obtained good agreement between the fit and the data 
(Figure 2). Goodness of fit quantification by the R 2 value 



(R 2 = \- 



-, where the jc, are the data points, y is 



the mean value of the data and the J[if)'s are the values of the 
model at times f,) gave an excellent score of R 2 = 0.99. 

Assuming that differences in growth between the primary tumor 
and its metastases should arise from interactions with the 
microenvironment, we fixed the proliferation parameter a for 
the metastases to the value obtained for the primary tumor 
growth. The only remaining parameters to be fixed were then b 
(driving the angiogenic stimulation) and m (controlling metastatic 
dissemination), allowing us to minimize the overall parameteriza- 
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Figure 2. Primary tumor growth. Comparison of the fit of the model 
and the data from [56]. Data are mean + standard error. 
doi:10.1371/journal.pone.0084249.g002 



tion of the model (two parameters for two data points). These last 
two coefficients were determined by fitting the model to the 
experimental metastatic data of Huang et al. [56] , the results of 
which are reported in Table 2. We obtained good agreement to 
the number and mean size of metastases. It was determined that 
with the estimated value of m a tumor of 200 mm' spreads a new 
metastasis every 0.77 days. 

The parameter estimation that we performed allowed us to 
simulate the experiment of [56] by using the parameters resulting 
from the model's fit (and /= 0). This gave further insight beyond 
the mere availability data could provide on the time development 
dynamics of the metastases and their final size distribution. 
Figure 3. A shows the growth in time of the total metastatic burden, 
while Figure 3.B depicts the colony size distribution at T = 32 days 
for an in silico replicate of the experiment performed in [56] . It 
reveals a nontrivial size distribution of the final metastatic colonies 
with a mode between 0.01 and 0.1 mm , and only one tumor with 
size larger than 1 0 mm' . At this time the total lung metastatic 
burden is 63.5 mm 3 distributed between 48.5 tumors. Simulation 
performed with a non-zero / and a value for/) extracted from [26] 
(see above) presented no significant difference in this setting 
compared to the simulation with /=0, hence justifying a posteriori 
our assumption of negligible effect of SIA in the setting of [56]. 

These simulations and parameter estimation show that our 
mathematical model is a possible theory describing growth and 
development of primary and secondary tumors in a 4T 1 cell line 
model. While in this context, systemic inhibition of angiogenesis had 
no significant importance in the small time range, it was concluded 
more broadly that the model, endowed with adequate parameter 
values, should be a useful theoretical tool for investigating a range of 
situations beyond the experimental setting of [56]. 

Simulation of the cancer history from the first cancer cell 
predicts uncontrolled metastatic burden 

Using our model and based on the parameters estimated in the 
previous section (Table 1), we were able to extrapolate to a totally 
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Table 2. Metastatic outputs. 







Value from [56] 


Computed by the model 


Median number of metastases (range) 


43 (4-107) 


43.03 


Mean size of metastases in mm 3 (range) 


1.47 (1.30-1.66) 


1.476 



Comparison of the fit of the model and the data from [56]. For the number of metastases, the reported model value is the number of tumors above a minimal visible 
size that we took to be 10 cells (tumors were counted using a dissecting microscope in [56]). Mean size was given as diameter in [56] and was converted here into 
volume using V= | xDx vv 2 , w= \ D. 
doi:1 0.1 371 /joumal.pone.0084249.t002 



new setting where the primary tumor starts with one cell instead of 
an already large number of cells (approximately 10 ). In so doing, 
we were able to simulate the whole cancer history, starting from 
one initial cancerous cell (and initial carrying capacity of 1 mm ) 
until the metastatic burden reached 5000 mm 3 , a burden we 
considered potentially lethal for a mouse. The simulation 
predicted this would happen 62.7 days after the first primary 
tumor cancer cell. Time development of the primary tumor 
volume, metastatic burden, total number and mean size of 
metastases as well as inhibitor amount in the host are plotted in 
Figure 4. 

Interestingly, the model simulation predicts that the metastatic 
burden would overcome the primary tumor mass, implying that 
the mouse would probably die from growth of its secondary lesions 
rather than from the initial tumor. This is consistent with the 
metastatically aggressive phenotype of the 4T1 cell line. Quanti- 
fication of the number of metastases reveals a final number of 
about 217 tumors, lots of them being small (Figure 4C) and 
probably undetectable in an experimental setting. Simulation with 
the same set of parameters but neglecting the effect of SIA (I— 0) 
showed no detectable difference on this time frame. Significant 
changes are observed later on, for volumes that are not considered 
to be physiologically relevant. This confirms that for the 4T1 cell 
line, metastases do develop and do not exhibit global dormancy, 
even when SIA is present with the inhibitor production parameter 
value extracted from [26]. Thus, based on biologically relevant 
parameters, our simulation results suggest large growth of the 



70 




Days 

A. Metastatic burden 

Figure 3. In silico simulation of the experiment from [56]. A. Time 
end time T = 32 days (log-scale on the x-axis). 
doi:1 0.1 371 /journal. pone.0084249.g003 



metastatic burden for the 4T 1 cell line when starting from the first 
cancer cell, with a final metastatic volume larger than the primary 
tumor. 

Higher production of systemic angiogenesis inhibitor 
could result in long-term stable global dormancy in a 
population of self-inhibiting metastases 

The previous simulations used parameter values derived from 
experimental data of a situation were metastases do develop and 
grow, because this is the only case where metastases are 
measurable and data are available. However we are interested 
here in global dormancy and situations where the metastatic 
population could remain ultimately small. We postulate that this 
could happen when production of the angiogenesis inhibitor, 
represented by parameter p in our model, is significantly higher. 
Simulation results plotted in Figure 5 were obtained using a value 
p = 2.5x10 4 mg'mm 3 day , i.e., a value about 30 times that 
extracted from [26]. From our previous modeling analysis and 
formula (3), higher production of inhibitor also proportionally 
increases the local inhibition parameters d and dp. In the 
simulation reported in Figure 5, we kept all the other parameter 
values unchanged from Table 1 and fixed the initial primary 
tumor volume to ^(0) = 1 cell and the initial primary tumor 
carrying capacity to Kp(0) — 1 mm . We simulated the system over 
a time of 350 days, covering the estimated lifespan of a mouse after 
appearance of an initial malignant cell. We focused on asymptotic 
behavior and possible convergence of the system to a steady state. 




Size 



B. Colonies size distribution at the final time 
levelopment of the metastatic burden. B. Colonies size distribution at the 
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Figure 4. Simulation of the cancer history from the first cancer cell. Parameter values are the ones resulting from the fit to the data of [56], 

reported in Table 1. 
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In Figure 5, the primary tumor volume, number and total burden 
of the metastases, the time evolution of the global inhibitor 
quantity and the size distribution of the metastases at the end time, 
are plotted. 

In this context, the first cancer cell initiates the disease by 
growing and generating a first pool of metastases, but the 
metastatic burden then quickly overshadows the growth of the 
primary lesion (Figure 5.A). The primary tumor reaches a small 
maximal size of 21.2 mm' at time 82.9 days (Figure 5.A) and then 
shrinks due to inhibition of angiogenesis provoked by the distant 
metastases. There is a slowdown and eventual stabilization of the 
metastatic burden, with a plateau value of about 2200 mm . The 
burden is composed of a large number of metastases (Figure 5.B), 
most of them being occult micro-metastases as can be seen in the 
final size distribution (Figure 5.C). This interesting feature of the 
model simulation could be an in silica replicate of the aforemen- 
tioned situations of cancer without disease [5] . In our model it 
translates into an asymptotical steady state for the metastatic 



burden while it is still composed of small lesions. The general 
dynamics of the metastatic burden results from the balance of two 
stimulating forces; growth and spread of new individuals, 
competing with systemic inhibition of angiogenesis. Stimulation 
depends on the parameters a and b, which capture the growth 
process, and m and a, which capture spreading. Inhibition depends 
on e and k, as well as on p, which controls the value of d. The 
present values of the parameters generated long-term stabilization 
of the mass. The size distribution of the population of secondary 
tumors at time T = 350 days is revealed to be non trivial, with 
different numbers in the various size ranges. By 350 days, all the 
metastases had volume lower than 10 mm . 

In sum, assuming substantial systemic inhibition of angiogenesis, 
we theoretically obtained an in silico replicate of a situation in 
which an important population of dormant micro-metastases 
inhibiting each others' growths is present, with a possibly non- 
lethal final total metastatic burden. This situation was seen to 
result when a 30-fold higher value for the inhibitor production 
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Figure 5. Large time simulation for large inhibitor production (p= 2.5 10 4 mg'mm 3 day \ d=2.16 mm 2 day '). The model 
predicts stabilization of the metastatic burden to a situation where the whole metastatic population is in a global dormancy state. 
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parameter p was used, compared to the case of growth of a breast 
cancer line 4T1 extracted from the literature [26], where 
unlimited expansion of the total metastatic burden was forecast. 

Discussion 

We propose an organism-scale model for the development of a 
primary tumor and a population of secondary tumors that takes 
into account systemic inhibiting interactions among tumors due to 
the release of a circulating angiogenesis inhibitor. The model 
proves to be able to describe in vivo data of primary tumor and 
metastatic development and allows inference of information not 
revealed by the experimental data, including the size density 
distribution of metastases and their total number. Endowed with 
biologically relevant parameter values, our model is a potentially 
vital tool for the theoretical study of metastatic dynamics. 

It was used here to investigate the whole cancer history from the 
first cancer cell, and predicted that for the metastatically aggressive 



4T1 cell line, metastases would grow unbounded for a physiolog- 
ically relevant set of parameter values. In this case, the total 
metastatic burden was found to become larger than the primary 
tumor mass and probably would be responsible for death of the 
animal. SIA effects were seen to be negligible in this context. A 
higher production rate of the inhibitor, by contrast, could 
theoretically make the primary tumor appearance and growth 
only a transient event, giving way to a distinct process of tumor 
development where, due to eventual self-inhibition of angiogenesis 
at the organism scale, global dormancy is imposed on the entire 
tumor/metastasis system, stabilizing the cancer disease. Our 
analysis shows that SIA could conceivably create such a situation, 
although it would require a very high value of the inhibitor 
production rate - some 30-fold the value extracted from [26] - 
which does not appear to be physiological. This suggests that SIA 
alone is probably not sufficient to induce spontaneous global 
dormancy and that other processes (such as immune effects) are 
probably significantly involved. For now, however, our conclusions 
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are limited by the current lack of data on systemic inhibition of 
tumor development. A study of interactions among multiple tumor 
implants in controlled immune contexts could shed more light. 

Meanwhile, these results inform the human situation by 
providing elements of explanation for the high prevalence of 
occult tumors found in autopsy studies. The results as well could 
inform the consequences of chronic antiangiogenic intervention 
[58] . As a case in point, progression of cancer disease in 
individuals having a low production rate of inhibitor might be 
forestalled, even outside an outright cure, by chronic external 
administration of supplementary inhibitory agents that could 
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